#   Program:  R
#   Task:     Plots
#   Project:  Markets 
#   Stata:    MAR09-DescGraph
#   Date:     07-01-2024


setwd("R:/Current Research/Markets-Pridemore/Work")

#################################
#         PACKAGES              #
#################################
library("foreign")
library("ggplot2")
library("dplyr")
library("gridExtra")
library("tidyverse")
library("patchwork")
source("panel.cor.R")

#################################
#     Loading Data              #
#################################
data<-read.dta("MAR09-DescGraph-old.dta")
inter.data<-read.csv("00-Market-InteractionData.csv")
avg.data<-read.dta("avg2.dta")
avg.data<-data.frame(na.omit(avg.data))
################################################################################
#############################################################
#       Setting Up Data for Pairs plot              #
#############################################################
pairs.plot.data  <- data %>% select(lrhom, fraser, infantmort, edu, unemp,  
                                    popdense, perurban,  sexratio)
pairs.plot.data<-data.frame(na.omit(pairs.plot.data))
colnames(pairs.plot.data)<-c("Ln Homicide", "Fraser", "Infant Mort", 
                             "Education Index", "Unemployment", 
                             "Population Density", "% Urban", "Sex Ratio")
pairs.plot.data.2  <- data %>% select(lrhom, fraser, infantmort, edu, unemp, 
                                      popdense, perurban,  sexratio, rhom,
                                      heritage, rol, KOFGI,  trust, notrust)
pairs.plot.data.2<-data.frame(na.omit(pairs.plot.data.2))
colnames(pairs.plot.data.2)<-c("Ln Homicide", "Fraser", "Infant Mort", 
                             "Education Index", "Unemployment", 
                             "Population Density", "% Urban", "Sex Ratio",
                             "Homicide Rates", "Heritage", "ROL", "Global", 
                             "Trust", "No Trust")
pairs.plot.data.3  <- data %>% select(lrhom, fraser, infantmort, edu, unemp, 
                                      popdense, perurban,  sexratio, rhom,
                                      heritage, rol)

pairs.plot.data.3<-data.frame(na.omit(pairs.plot.data.3))
colnames(pairs.plot.data.3)<-c("Ln Homicide", "Fraser", "Infant Mort", 
                               "Education Index", "Unemployment", 
                               "Population Density", "% Urban", "Sex Ratio",
                               "Homicide Rates", "Heritage", "ROL")
#############################################################
#                           Pairs Plots             #
#############################################################
pairs(pairs.plot.data,  upper.panel = panel.smooth, lower.panel = panel.cor)
pairs(pairs.plot.data.2,  upper.panel = panel.smooth, lower.panel = panel.cor)
pairs(pairs.plot.data.3,  upper.panel = panel.smooth, lower.panel = panel.cor)

###############################################################################
#############################################################
#                Avg Plots            #
#############################################################
ggplot(avg.data, aes(x = year)) +
  geom_line(aes(y = rhom, color = "Homicide Rates", linetype = "Homicide Rates")) +
  geom_line(aes(y = fraser, color = "Fraser", linetype = "Fraser")) +
  scale_color_manual(name = "Legend",
                     values = c("Homicide Rates" = "black", "Fraser" = "black")) +  
  scale_linetype_manual(name = "Legend",
                        values = c("Homicide Rates" = "solid", "Fraser" = "dashed")) +
  labs(x = "Year", y = "Average Yearly Homicide Rates per 100k pop & Fraser") +
  scale_x_continuous(breaks = seq(2000, 2019, by = 1)) +
  theme_bw()


ggplot(avg.data, aes(x = year)) +
  geom_line(aes(y = rhom, color = "Homicide Rates", linetype = "Homicide Rates")) +
  geom_line(aes(y = fraser, color = "Fraser", linetype = "Fraser")) +
  scale_color_manual(name = "Legend",
                     values = c("Homicide Rates" = "black", "Fraser" = "black")) +  
  scale_linetype_manual(name = "Legend",
                        values = c("Homicide Rates" = "solid", "Fraser" = "dashed")) +
  labs(x = "Year", y = "Average Yearly Homicide Rates per 100k pop & Fraser") +
  scale_x_continuous(breaks = seq(2000, 2019, by = 1)) +
  theme_bw()+ 
  theme(panel.grid.major = element_blank(),  
        panel.grid.minor = element_blank())  

#########################################
#     All in one graph                  #
#########################################
ggplot(avg.data, aes(x = year)) +
  # Full sample
  geom_line(aes(y = rhom, color = "Homicide Rate (All)", linetype = "Homicide Rate (All)")) +
  geom_line(aes(y = fraser, color = "Fraser (All)", linetype = "Fraser (All)")) +
  
  # Below-median Fraser countries
  geom_line(aes(y = Lrhom, color = "Homicide ≤Median Fraser", linetype = "Homicide ≤Median Fraser")) +
  geom_line(aes(y = Lfraser, color = "Fraser ≤Median", linetype = "Fraser ≤Median")) +
  
  # Above-median Fraser countries
  geom_line(aes(y = Hrhom, color = "Homicide >Median Fraser", linetype = "Homicide >Median Fraser")) +
  geom_line(aes(y = Hfraser, color = "Fraser >Median", linetype = "Fraser >Median")) +
  
  # Manual greyscale styling
  scale_color_manual(name = "Legend",
                     values = c("Homicide Rate (All)" = "black",
                                "Fraser (All)" = "grey50",
                                "Homicide ≤Median Fraser" = "black",
                                "Fraser ≤Median" = "gray50",
                                "Homicide >Median Fraser" = "black",
                                "Fraser >Median" = "gray50")) +
  
  scale_linetype_manual(name = "Legend",
                        values = c("Homicide Rate (All)" = "solid",
                                   "Fraser (All)" = "dashed",
                                   "Homicide ≤Median Fraser" = "dotted",
                                   "Fraser ≤Median" = "dotdash",
                                   "Homicide >Median Fraser" = "longdash",
                                   "Fraser >Median" = "twodash")) +
  
  labs(x = "Year", y = "Homicide Rate (per 100k) & Fraser Index") +
  scale_x_continuous(breaks = seq(2000, 2019, by = 1)) +
  scale_y_continuous(breaks = seq(0, 12, by = 1)) +
  
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom")

############################################
#  3 seprate but comparable graphs        #
###########################################
# Plot 1: Full sample
p1 <- ggplot(avg.data, aes(x = year)) +
  geom_line(aes(y = rhom, linetype = "Homicide Rate (All)")) +
  geom_line(aes(y = fraser, linetype = "Fraser (All)")) +
  scale_linetype_manual(name = "All Nations",
                        values = c("Homicide Rate (All)" = "solid",
                                   "Fraser (All)" = "dashed")) +
  labs(title = "All Nations",
       x = "Year", y = "Homicide Rate & Fraser Index") +
  scale_y_continuous(breaks = seq(0, 12, 1), limits = c(0, 12))+
theme_bw() +
  theme(panel.grid = element_blank(), legend.position = "bottom")

# Plot 2: ≤Median Fraser countries
p2 <- ggplot(avg.data, aes(x = year)) +
  geom_line(aes(y = Lrhom, linetype = "Homicide ≤Median Fraser")) +
  geom_line(aes(y = Lfraser, linetype = "Fraser ≤Median")) +
  scale_linetype_manual(name = "Low Fraser",
                        values = c("Homicide ≤Median Fraser" = "solid",
                                   "Fraser ≤Median" = "dashed")) +
  labs(title = "Low Fraser (≤Median Fraser) Nations",
       x = "Year", y = "Homicide Rate & Fraser Index") +
  scale_y_continuous(breaks = seq(0, 12, 1), limits = c(0, 12))+
theme_bw() +
  theme(panel.grid = element_blank(), legend.position = "bottom")

# Plot 3: >Median Fraser countries
p3 <- ggplot(avg.data, aes(x = year)) +
  geom_line(aes(y = Hrhom, linetype = "Homicide >Median Fraser")) +
  geom_line(aes(y = Hfraser, linetype = "Fraser >Median")) +
  scale_linetype_manual(name = "High Fraser",
                        values = c("Homicide >Median Fraser" = "solid",
                                   "Fraser >Median" = "dashed")) +
  labs(title = "High (>Median Fraser) Nations",
       x = "Year", y = "Homicide Rate & Fraser Index") +
  scale_y_continuous(breaks = seq(0, 12, 1), limits = c(0, 12))+
theme_bw() +
  theme(panel.grid = element_blank(), legend.position = "bottom")

combined_plot <- (p1 + p2 + p3 ) +
  patchwork::plot_layout(ncol = 2)
  
# Print or export
print(combined_plot)


###############################################################################
#################################################################
#                   Predicted Plots  Set up                    #
#################################################################
fraser.data <- inter.data %>% 
  select(fraser, yhat_fraserdir, yhat_lowfraser, yhat_highfraser) %>%
  na.omit() %>%
  rename(
    Direct = yhat_fraserdir,
    Low = yhat_lowfraser,
    High = yhat_highfraser
  )

#################################################################
#                   Predicted Plots  Fraser, Low, and High      #
#################################################################
ggplot(fraser.data) +
  geom_line(aes(x = fraser, y = Direct, color = "Direct"), linetype = "solid") +
  geom_line(aes(x = fraser, y = Low, color = "Low"), linetype = "dashed") +
  geom_line(aes(x = fraser, y = High, color = "High"), linetype = "dotted") +
  labs(x = "Fraser", y = "Predicted Ln Homicide Rates per 100k pop") +
  scale_x_continuous(name = "Fraser") +
  scale_y_continuous(name = "Predicted Ln Homicide Rates per 100k pop") +
  scale_color_manual(name = "Legend",  # Legend title
                     labels = c("Fraser", "Low", "High"),  # Legend labels
                     values = c("black", "black", "black")) +  # Colors for lines
  theme_bw()

ggplot(fraser.data) +
  geom_line(aes(x = fraser, y = Direct, color = "Direct"), linetype = "solid") +
  geom_line(aes(x = fraser, y = Low, color = "Low"), linetype = "dashed") +
  geom_line(aes(x = fraser, y = High, color = "High"), linetype = "dotted") +
  labs(x = "Fraser", y = "Predicted Ln Homicide Rates per 100k pop") +
  scale_x_continuous(name = "Fraser") +
  scale_y_continuous(name = "Predicted Ln Homicide Rates per 100k pop") +
  scale_color_manual(name = "Legend",  # Legend title
                     labels = c("Fraser", "Low", "High"),  # Legend labels
                     values = c("black", "black", "black")) +  # Colors for lines
  theme_bw()+
  theme(panel.grid.major = element_blank(),  
        panel.grid.minor = element_blank())  
################################################################################
#######################################################
#             INTERACTION PLOT SET UP                 #
#######################################################
colnames(inter.data)<-c("fraser", "infantmort", "yhat.infant", "lgdp", 
                        "yhat.gdp.fraser","yhat.gdp", "yhat.fras.dir",
                        "yhat.low.fraser", "yhat.high.fraser")
interaction<-inter.data  %>% 
  select(fraser, infantmort, yhat.infant)
interaction.2 <- subset(inter.data[inter.data$fraser == 2, c("fraser", "infantmort", "yhat.infant")])
interaction.8 <- subset(inter.data[inter.data$fraser == 8, c("fraser", "infantmort", "yhat.infant")])
interaction_combined <- rbind(transform(interaction.2, Group = "Low"), 
                              transform(interaction.8, Group = "High"))

gdp.inter.data <- tibble(
  fraser = c(2, 2, 2, 2, 8, 8, 8, 8),
  lgdp = c(6, 8, 10, 12, 6, 8, 10, 12),
  yhat = c(10.232481, 5.0527908, -0.12689962, -5.3065901, 1.8578305, 1.2131263, 0.56842211, -0.0762821),
  Group = c("Low", "Low", "Low", "Low", "High", "High", "High", "High")
)

#########################################################
#               First Plot InfantX Fraser               #
#########################################################
infant.plot <- ggplot(data = interaction_combined, aes(x = infantmort, y = yhat.infant, color = Group, linetype = Group)) +
  geom_line() +
  labs(x = "Infant Mortality", y = "Predicted Ln Homicide Rates per 100k pop", title = "Infant Mortality X Fraser") +
  scale_color_manual(name = "Fraser",  # Change legend title to Fraser
                     values = c("black", "black"), 
                     labels = c("Low", "High")) +  # Adjust labels as per your data
  scale_linetype_manual(name = "Fraser",  # Change legend title to Fraser
                        values = c("solid", "dashed"), 
                        labels = c("Low", "High")) +  # Adjust labels as per your data
  ylim(range(c(interaction_combined$yhat.infant, gdp.inter.data$yhat))) +  # Set y-axis limits based on combined data range
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, margin = margin(b = 20)))  # Center align title

infant.plot.n <- ggplot(data = interaction_combined, aes(x = infantmort, y = yhat.infant, color = Group, linetype = Group)) +
  geom_line() +
  labs(x = "Infant Mortality", y = "Predicted Ln Homicide Rates per 100k pop", title = "Infant Mortality X Fraser") +
  scale_color_manual(name = "Fraser",  # Change legend title to Fraser
                     values = c("black", "black"), 
                     labels = c("Low", "High")) +  # Adjust labels as per your data
  scale_linetype_manual(name = "Fraser",  # Change legend title to Fraser
                        values = c("solid", "dashed"), 
                        labels = c("Low", "High")) +  # Adjust labels as per your data
  ylim(range(c(interaction_combined$yhat.infant, gdp.inter.data$yhat))) +  # Set y-axis limits based on combined data range
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, margin = margin(b = 20))) +
  theme(panel.grid.major = element_blank(),  
        panel.grid.minor = element_blank())  




#########################################################
#               Second Plot LN GDP X Fraser               #
#########################################################
gdp.plot <- ggplot(data = gdp.inter.data, aes(x = lgdp, y = yhat, color = Group, linetype = Group)) +
  geom_line() +
  labs(x = "Ln GDP", y = "Predicted Ln Homicide Rates per 100k pop", title = "Ln GDP X Fraser") +
  scale_color_manual(name = "Fraser",  # Change legend title to Fraser
                     values = c("black", "black"), 
                     labels = c("Low", "High")) +  # Adjust labels as per your data
  scale_linetype_manual(name = "Fraser",  # Change legend title to Fraser
                        values = c("solid", "dashed"), 
                        labels = c("Low", "High")) +  # Adjust labels as per your data
  ylim(range(c(interaction_combined$yhat.infant, gdp.inter.data$yhat))) +  # Set y-axis limits based on combined data range
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, margin = margin(b = 20)))  # Center align title

gdp.plot.n <- ggplot(data = gdp.inter.data, aes(x = lgdp, y = yhat, color = Group, linetype = Group)) +
  geom_line() +
  labs(x = "Ln GDP", y = "Predicted Ln Homicide Rates per 100k pop", title = "Ln GDP X Fraser") +
  scale_color_manual(name = "Fraser",  # Change legend title to Fraser
                     values = c("black", "black"), 
                     labels = c("Low", "High")) +  # Adjust labels as per your data
  scale_linetype_manual(name = "Fraser",  # Change legend title to Fraser
                        values = c("solid", "dashed"), 
                        labels = c("Low", "High")) +  # Adjust labels as per your data
  ylim(range(c(interaction_combined$yhat.infant, gdp.inter.data$yhat))) +  # Set y-axis limits based on combined data range
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, margin = margin(b = 20))) +
  theme(panel.grid.major = element_blank(),  
        panel.grid.minor = element_blank())  




##############################################
# Combine plots with shared legend
###########################################
combined_plot <- infant.plot + gdp.plot +
  plot_layout(guides = 'collect')  # Ensures legends are combined

# Display the combined plot
combined_plot


combined_plot.n <- infant.plot.n + gdp.plot.n +
  plot_layout(guides = 'collect')  # Ensures legends are combined
combined_plot.n


################################################################################
#                 CURVLINEAR                                          #
###############################################################################
curv<-read.dta("ppFraser2.dta")
ggplot(curv, aes(x = at1, y = margin)) +
  geom_line(color = "black", linetype = "solid") +
  geom_ribbon(aes(ymin = ci_lb, ymax = ci_ub), alpha = 0.2, fill = "gray50") +
  labs(
    x = "Fraser Index",
    y = "Predicted Ln Homicide Rate",
    title = " "   ) +
  scale_x_continuous(breaks = seq(min(curv$at1), max(curv$at1), 1)) +
  scale_y_continuous(breaks = seq(0, 10, 1)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())


ggplot(curv, aes(x = at1, y = margin)) +
  geom_line(color = "black", linetype = "solid") +
  geom_ribbon(aes(ymin = ci_lb, ymax = ci_ub), alpha = 0.2, fill = "gray50") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(
    x = "Fraser Index",
    y = "Predicted Ln Homicide Rate",
    title = " "
  ) +
  scale_x_continuous(breaks = seq(min(curv$at1), max(curv$at1), 1)) +
  scale_y_continuous(breaks = seq(0, 10, 1)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())
